# Author: Gregory Smith
# Date: 12 - 3 - 2018
# Title: Secret But Constrained Appendix B
# Job: This File Includes all of the Robustness Checks Included in the Final
# Appendix that rely upon bootstrapped clustered SEs 

# Note: The script smith-2019-covert-final-tidy-output-functions needs to be run
#  for this code to work

# Clear the Workspace Before Running 

#This file takes a significant amount of time to run. Paramater estimates
# and goodness of fit statistics have also been saved as tidy data frames. 

#---- Packages 

#---- Packages 

list.of.packages <- c("tidyverse")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

library(tidyverse)

# install_version("clusterSEs", version = "2.5", repos = "http://cran.us.r-project.org")
library(clusterSEs)  # Clustered Standard Errors 

#---- Set the Seed 

set.seed(175)

#------ Loads the Data 

dta <- read_csv("data/smith_2019_secret_but_constrained_dta.csv")

#------ Pooled MIDs:

# Note: the following code assesses the effects of divided government on 
# pooled MIDs and re-runs models with major and minor MIDs as the dv with an 
# alternatve peace year specificaiton that drops peace years when either


pooled_mid_base <- glm(us_mid_onset ~  divided + 
                         factor(president) - 1 +
                         us_onset_py + us_onset_py_2  + us_onset_py_3,
                       data = dta, family = "binomial")

pooled_mid_base_cse <- cluster.bs.glm(pooled_mid_base, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(pooled_mid_base_cse, pooled_mid_base, model_name = "pooled_mid_base")

pooled_mid_full <- glm(us_mid_onset ~  divided + 
                         election  + approval_q1 +
                         democracy + cinc_proportion + clean_russ_inf + 
                         oversight_change + 
                         factor(president) - 1 +
                         us_onset_py + us_onset_py_2  + us_onset_py_3,
                       data = dta, family = "binomial")

pooled_mid_cse <- cluster.bs.glm(pooled_mid_full, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(pooled_mid_cse, pooled_mid_full, model_name = "pooled_mid_full")


#------- Major and Minor MIDs with Years Dropped When Major/Minor Ongoing 

minor_mid_test_base <- glm(minor_us_onset ~  divided + 
                             factor(president) - 1 +
                             us_onset_py + us_onset_py_2  + us_onset_py_3,
                           data = dta, family = "binomial")


minor_mid_test_base_cse <- cluster.bs.glm(minor_mid_test_base, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_mid_test_base_cse, minor_mid_test_base, model_name = "minor_mid_test_base")


minor_mid_test_full <- glm(minor_us_onset ~  divided + 
                             election  + approval_q1 +
                             democracy + cinc_proportion + clean_russ_inf + 
                             oversight_change + 
                             factor(president) - 1 +
                             us_onset_py + us_onset_py_2  + us_onset_py_3,
                           data = dta, family = "binomial")

minor_mid_test_full_cse <- cluster.bs.glm(minor_mid_test_full, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_mid_test_full_cse, minor_mid_test_full, model_name = "minor_mid_test_full")


major_mid_test_base <- glm(major_us_onset ~  divided + 
                             factor(president) - 1 +
                             us_onset_py + us_onset_py_2  + us_onset_py_3,
                           data = dta, family = "binomial")

major_mid_test_base_cse <- cluster.bs.glm(major_mid_test_base, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_mid_test_base_cse, major_mid_test_base, model_name = "major_mid_test_base")


major_mid_test_full <- glm(major_us_onset ~  divided + 
                             election  + approval_q1 +
                             democracy + cinc_proportion + clean_russ_inf + 
                             oversight_change + 
                             factor(president) - 1 +
                             us_onset_py + us_onset_py_2  + us_onset_py_3,
                           data = dta, family = "binomial")

major_mid_test_full_cse <- cluster.bs.glm(major_mid_test_full, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_mid_test_full_cse, major_mid_test_full, model_name = "major_mid_test_full")

#----- Alternative Measures of Party Control (Majority Pres. Party and LPPC Scores)

# O'Rourke Covert Operations 

majority_orke_base <- glm(orourke_onset ~  prezmajority + 
                            factor(president) - 1  +
                            orourke_peace_years + orourke_peace_years_2  + orourke_peace_years_3,
                          data = dta, family = "binomial")

majority_orke_base_cse <- cluster.bs.glm(majority_orke_base, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(majority_orke_base_cse, majority_orke_base, model_name = "majority_orke_base")


majority_orke_full <- glm(orourke_onset ~  prezmajority + 
                            election  + approval_q1 +
                            democracy + cinc_proportion + clean_russ_inf + 
                            oversight_change + 
                            factor(president) - 1 +
                            orourke_peace_years + orourke_peace_years_2  + orourke_peace_years_3,
                          data = dta, family = "binomial")

majority_orke_full_cse <- cluster.bs.glm(majority_orke_full, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(majority_orke_full_cse, majority_orke_full, model_name = "majority_orke_full")


lppc_orke_base <- glm(orourke_onset ~  lppc + 
                        factor(president) - 1 +
                        orourke_peace_years + orourke_peace_years_2  + orourke_peace_years_3,
                      data = dta, family = "binomial")

lppc_orke_base_cse <- cluster.bs.glm(lppc_orke_base, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(lppc_orke_base_cse, lppc_orke_base, model_name = "lppc_orke_base")


lppc_orke_full <- glm(orourke_onset ~  lppc + 
                        election  + approval_q1 +
                        democracy + cinc_proportion + clean_russ_inf + 
                        oversight_change + 
                        factor(president) - 1 +
                        orourke_peace_years + orourke_peace_years_2  + orourke_peace_years_3,
                      data = dta, family = "binomial")

lppc_orke_full_cse <- cluster.bs.glm(lppc_orke_full, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(lppc_orke_full_cse, lppc_orke_full, model_name = "lppc_orke_full")


#---- Models that use the Pooled (O'Rourke + Berger) Sample 

majority_cmb_base <- glm(us_covert_onset ~  prezmajority + 
                           factor(president) - 1 +
                           us_peace_years + us_peace_years_2 + us_peace_years_3, 
                         data = dta, family = "binomial")

majority_cmb_base_cse <- cluster.bs.glm(majority_cmb_base, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(majority_cmb_base_cse, majority_cmb_base, model_name = "majority_cmb_base")


majority_cmb_full <- glm(us_covert_onset ~  prezmajority + 
                           election  + approval_q1 +
                           democracy + cinc_proportion + clean_russ_inf + 
                           oversight_change +  
                           factor(president) - 1 +
                           us_peace_years + us_peace_years_2 + us_peace_years_3, 
                         data = dta, family = "binomial")

majority_cmb_full_cse <- cluster.bs.glm(majority_cmb_full, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(majority_cmb_full_cse, majority_cmb_full, model_name = "majority_cmb_full")


lppc_cmb_base <- glm(us_covert_onset ~  lppc + 
                       factor(president) - 1 +
                       us_peace_years + us_peace_years_2 + us_peace_years_3, 
                     data = dta, family = "binomial")

lppc_cmb_base_cse <- cluster.bs.glm(lppc_cmb_base, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(lppc_cmb_base_cse, lppc_cmb_base, model_name = "lppc_cmb_base")


lppc_cmb_full <- glm(us_covert_onset ~  lppc + 
                       election  + approval_q1 +
                       democracy + cinc_proportion + clean_russ_inf + 
                       oversight_change +  
                       factor(president) - 1 +
                       us_peace_years + us_peace_years_2 + us_peace_years_3, 
                     data = dta, family = "binomial")

lppc_cmb_full_cse <- cluster.bs.glm(lppc_cmb_full, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(lppc_cmb_full_cse, lppc_cmb_full, model_name = "lppc_cmb_full")



#----- Minor US Initiated GML MIDs 

minor_us_base_mid_mjt <- glm(minor_us_onset ~  prezmajority + 
                               factor(president) - 1 + 
                               minor_us_py + I(minor_us_py^2)  +
                               I(minor_us_py^3), 
                             data = dta, family = "binomial")

minor_us_base_mid_mjt_cse <- cluster.bs.glm(minor_us_base_mid_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_base_mid_mjt_cse, minor_us_base_mid_mjt, model_name = "minor_us_base_mid_mjt")


minor_us_mid_mjt <- glm(minor_us_onset ~  prezmajority + 
                          election  + approval_q1 +
                          democracy + cinc_proportion + clean_russ_inf + 
                          oversight_change + 
                          factor(president) - 1 + 
                          minor_us_py + I(minor_us_py^2)  +
                          I(minor_us_py^3), 
                        data = dta, family = "binomial")

minor_us_mid_mjt_cse <- cluster.bs.glm(minor_us_mid_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_mid_mjt_cse, minor_us_mid_mjt, model_name = "minor_us_mid_mjt")


minor_us_base_mid_lppc <- glm(minor_us_onset ~  lppc + 
                                factor(president) - 1 + 
                                minor_us_py + I(minor_us_py^2)  +
                                I(minor_us_py^3), 
                              data = dta, family = "binomial")

minor_us_base_mid_lppc_cse <- cluster.bs.glm(minor_us_base_mid_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_base_mid_lppc_cse, minor_us_base_mid_lppc, model_name = "minor_us_base_mid_lppc")


minor_us_mid_lppc <- glm(minor_us_onset ~  lppc + 
                           election  + approval_q1 +
                           democracy + cinc_proportion + clean_russ_inf + 
                           oversight_change + 
                           factor(president) - 1 + 
                           minor_us_py + I(minor_us_py^2)  +
                           I(minor_us_py^3), 
                         data = dta, family = "binomial")

minor_us_mid_lppc_cse <- cluster.bs.glm(minor_us_mid_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_mid_lppc_cse, minor_us_mid_lppc, model_name = "minor_us_mid_lppc")

#----- Major US Initiated GML MIDs 

major_us_base_mid_mjt <- glm(major_us_onset ~  prezmajority + 
                               factor(president) - 1 + 
                               major_us_py + I(major_us_py^2)  +
                               I(major_us_py^3), 
                             data = dta, family = "binomial")

major_us_base_mid_mjt_cse <- cluster.bs.glm(major_us_base_mid_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_base_mid_mjt_cse, major_us_base_mid_mjt, model_name = "major_us_base_mid_mjt")


major_us_mid_mjt <- glm(major_us_onset ~  prezmajority + 
                          election  + approval_q1 +
                          democracy + cinc_proportion + clean_russ_inf + 
                          oversight_change + 
                          factor(president) - 1 + 
                          major_us_py + I(major_us_py^2)  +
                          I(major_us_py^3), 
                        data = dta, family = "binomial")

major_us_mid_mjt_cse <- cluster.bs.glm(major_us_mid_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_mid_mjt_cse, major_us_mid_mjt, model_name = "major_us_mid_mjt")


major_us_base_mid_lppc <- glm(major_us_onset ~  lppc + 
                                factor(president) - 1 + 
                                major_us_py + I(major_us_py^2)  +
                                I(major_us_py^3), 
                              data = dta, family = "binomial")

major_us_base_mid_lppc_cse <- cluster.bs.glm(major_us_base_mid_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_base_mid_lppc_cse, major_us_base_mid_lppc, model_name = "major_us_base_mid_lppc")


major_us_mid_lppc <- glm(major_us_onset ~  lppc + 
                           election  + approval_q1 +
                           democracy + cinc_proportion + clean_russ_inf + 
                           oversight_change + 
                           factor(president) - 1 + 
                           major_us_py + I(major_us_py^2)  +
                           I(major_us_py^3), 
                         data = dta, family = "binomial")

major_us_mid_lppc_cse <- cluster.bs.glm(major_us_mid_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_mid_lppc_cse, major_us_mid_lppc, model_name = "major_us_mid_lppc")


#----- The Following Models Rely Upon Fordham and Sarver's Measure of Force 

#  Minor Uses of Force 

minor_us_base_force_div <- glm(minor_force_onset ~  divided + 
                                 factor(president) - 1 + 
                                 minor_force_peace + I(minor_force_peace^2)  +
                                 I(minor_force_peace^3), 
                               data = dta, family = "binomial")

minor_us_base_force_div_cse <- cluster.bs.glm(minor_us_base_force_div, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_base_force_div_cse, minor_us_base_force_div, model_name = "minor_us_base_force_div")


minor_us_force_div <- glm(minor_force_onset ~  divided + 
                            election  + approval_q1 + 
                            democracy + cinc_proportion + clean_russ_inf + 
                            oversight_change + 
                            factor(president) - 1 + 
                            minor_force_peace + I(minor_force_peace^2)  +
                            I(minor_force_peace^3), 
                          data = dta, family = "binomial")

minor_us_force_div_cse <- cluster.bs.glm(minor_us_force_div, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_force_div_cse, minor_us_force_div, model_name = "minor_us_force_div")


minor_us_base_force_mjt <- glm(minor_force_onset ~  prezmajority + 
                                 factor(president) - 1 + 
                                 minor_force_peace + I(minor_force_peace^2)  +
                                 I(minor_force_peace^3), 
                               data = dta, family = "binomial")

minor_us_base_force_mjt_cse <- cluster.bs.glm(minor_us_base_force_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_base_force_mjt_cse, minor_us_base_force_mjt, model_name = "minor_us_base_force_mjt")


minor_us_force_mjt <- glm(minor_force_onset ~  prezmajority + 
                            election  + approval_q1 +
                            democracy + cinc_proportion + clean_russ_inf + 
                            oversight_change + 
                            factor(president) - 1 + 
                            minor_force_peace + I(minor_force_peace^2)  +
                            I(minor_force_peace^3), 
                          data = dta, family = "binomial")

minor_us_force_mjt_cse <- cluster.bs.glm(minor_us_force_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_force_mjt_cse, minor_us_force_mjt, model_name = "minor_us_force_mjt")


minor_us_base_force_lppc <- glm(minor_force_onset ~  lppc + 
                                  factor(president) - 1 + 
                                  minor_force_peace + I(minor_force_peace^2)  +
                                  I(minor_force_peace^3), 
                                data = dta, family = "binomial")

minor_us_base_force_lppc_cse <- cluster.bs.glm(minor_us_base_force_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_base_force_lppc_cse, minor_us_base_force_lppc, model_name = "minor_us_base_force_lppc")


minor_us_force_lppc <- glm(minor_force_onset ~  lppc + 
                             election  + approval_q1 +
                             democracy + cinc_proportion + clean_russ_inf + 
                             oversight_change + 
                             factor(president) - 1 + 
                             minor_force_peace + I(minor_force_peace^2)  +
                             I(minor_force_peace^3), 
                           data = dta, family = "binomial")

minor_us_force_lppc_cse <- cluster.bs.glm(minor_us_force_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(minor_us_force_lppc_cse, minor_us_force_lppc, model_name = "minor_us_force_lppc")

# Major Uses of Force 

major_us_base_force_div <- glm(major_force_onset ~  divided + 
                                 factor(president) - 1 + 
                                 major_force_peace + I(major_force_peace^2)  +
                                 I(major_force_peace^3), 
                               data = dta, family = "binomial")

major_us_base_force_div_cse <- cluster.bs.glm(major_us_base_force_div, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_base_force_div_cse, major_us_base_force_div, model_name = "major_us_base_force_div")


major_us_force_div <- glm(major_force_onset ~  divided + 
                            election  + approval_q1 +
                            democracy + cinc_proportion + clean_russ_inf + 
                            oversight_change + 
                            factor(president) - 1 + 
                            major_force_peace + I(major_force_peace^2)  +
                            I(major_force_peace^3), 
                          data = dta, family = "binomial")

major_us_force_div_cse <- cluster.bs.glm(major_us_force_div, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_force_div_cse, major_us_force_div, model_name = "major_us_force_div")


major_us_base_force_mjt <- glm(major_force_onset ~  prezmajority + 
                                 factor(president) - 1 + 
                                 major_force_peace + I(major_force_peace^2)  +
                                 I(major_force_peace^3), 
                               data = dta, family = "binomial")

major_us_base_force_mjt_cse <- cluster.bs.glm(major_us_base_force_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_base_force_mjt_cse, major_us_base_force_mjt, model_name = "major_us_base_force_mjt")


major_us_force_mjt <- glm(major_force_onset ~  prezmajority + 
                            election  + approval_q1 +
                            democracy + cinc_proportion + clean_russ_inf + 
                            oversight_change + 
                            factor(president) - 1 + 
                            major_force_peace + I(major_force_peace^2)  +
                            I(major_force_peace^3), 
                          data = dta, family = "binomial")

major_us_force_mjt_cse <- cluster.bs.glm(major_us_force_mjt, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_force_mjt_cse, major_us_force_mjt, model_name = "major_us_force_mjt")


major_us_base_force_lppc <- glm(major_force_onset ~  lppc + 
                                  factor(president) - 1 + 
                                  major_force_peace + I(major_force_peace^2)  +
                                  I(major_force_peace^3), 
                                data = dta, family = "binomial")

major_us_base_force_lppc_cse <- cluster.bs.glm(major_us_base_force_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_base_force_lppc_cse, major_us_base_force_lppc, model_name = "major_us_base_force_lppc")


major_us_force_lppc <- glm(major_force_onset ~  lppc + 
                             election  + approval_q1 +
                             democracy + cinc_proportion + clean_russ_inf + 
                             oversight_change + 
                             factor(president) - 1 + 
                             major_force_peace + I(major_force_peace^2)  +
                             I(major_force_peace^3), 
                           data = dta, family = "binomial")

major_us_force_lppc_cse <- cluster.bs.glm(major_us_force_lppc, dta, ~ ccode_current, boot.reps = 2000)
clstr_se_clean(major_us_force_lppc_cse, major_us_force_lppc, model_name = "major_us_force_lppc")

